%gridding with maximum between-class variance method(Otsu)
%p:input image , e.g. GSM16101_CH1-1.tif
% for a series of image processing, the "Drawgrid" had been noted,if you want to see the gridding result
%you can delete the note "%"
function para=otsugrid(p)
para=[];
st=cputime;
%read image
A=imread(p);
A=double(A);
%transfer 16 bit to 8 bit
A=A./256;
[h,w,d]=size(A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%horizontal projection
for i=1:h
    m=0;
    for j=1:w
          m=m+A(i,j);
    end
    Hy(i)=m/w;
end
%mean
hym=sum(Hy)/h;
%reduce mean
Hy1=Hy-hym;
%morphological reconstruction
Hy2=imreconstruct(Hy1,Hy);
clear Hy1;
%get the residue signal
Hy3=Hy-Hy2;
clear Hy2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hy3=uint8(Hy3);
%get threshold th by using the Otsu method
th=thresh_md(Hy3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%vertical projection
for i=1:w
    n=0;
    for j=1:h
          n=n+A(j,i);
    end
    Vx(i)=n/h;
end
%mean
vxm=sum(Vx)/w;
%reduce mean
Vx1=Vx-vxm;
%morphological reconstruction
Vx2=imreconstruct(Vx1,Vx);
clear Vx1;
%get the residue signal
Vx3=Vx-Vx2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vx3=uint8(Vx3);
%figure,plot(1:w,Vx);
%get the threshold th1 by using Otsu method
th1=thresh_md(Vx3);
th
th1
th2=uint8((th+th1)/2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%transfer the projection signal to binary signal
for i=1:h 
    if Hy3(i)>th 
          Hy1(i)=0; 
    else 
          Hy1(i)=255; 
    end 
end
for i=1:w 
    if Vx3(i)>th1
          Vx1(i)=0; 
    else 
          Vx1(i)=255; 
    end 
end
%figure ,plot(1:h,Hy1);
%figure,plot(1:w,Vx1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find coordinate y
hf=0;%flag
flag=0;
%find the region of 1 or 255
for i=1:h-2 %find the original pixel 
    if flag==0
        if Hy1(i)==255 && Hy1(i+1)==255 && Hy1(i+2)==255
         hf=hf+1;
         flag=1;
         hspotcoff(hf)=i;
         t=0;
       end
    end
    if flag==1%find the terminal pixel
        t=t+1;
       if Hy1(i+1)==0 &&t>4
           hf=hf+1;
           hspotcoff(hf)=i;
           flag=0;
       end
    end
end
clear Hy1;
%compute the middle of two grid lines
disp('compute the middle to get final horizontal grid lines');
hspotcoff1(1)=hspotcoff(2);
for i=3:2:hf-1
    hspotcoff1((i+1)/2)=round((hspotcoff(i)+hspotcoff(i+1))/2);
end
hspotcoff1(round((hf-1)/2+1))=hspotcoff(hf);
%figure ,plot(1:w,Vx1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%find cooridinate x
vf=0;%flag
flag=0;
for i=1:w-2
    if flag==0
        if Vx1(i)==255 &&Vx1(i+1)==255 && Vx1(i+2)==255
         vf=vf+1;
         flag=1;
        vspotcoff(vf)=i;
        t=0;
       end
    end
    if flag==1
        t=t+1;
       if Vx1(i+1)==0 &&t>4
           vf=vf+1;
           vspotcoff(vf)=i;
           flag=0;
       end
    end
end
clear Vx1;
%compute the middle of two grid lines
disp('compute the middle to get final vertical grid lines');
vspotcoff1(1)=vspotcoff(2);
for i=3:2:vf-1
    vspotcoff1((i+1)/2)=round((vspotcoff(i)+vspotcoff(i+1))/2);
end
vspotcoff1(round((vf-1)/2+1))=vspotcoff(vf);
%Drawgrid(p,h,w,hspotcoff,vspotcoff,hf,vf);
%reduce the first grid line which was for location
hf=hf-1;
vf=vf-1;
clear hspotcoff;
clear vspotcoff;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
et=cputime-st;
disp(['time for gridding',num2str(et)]);
%draw the original grid line
[s,hf1]=size(hspotcoff1);
[s,vf1]=size(vspotcoff1);
mm=5;
hspotcoff1(1)=hspotcoff1(1)-mm;
if hspotcoff1(1)<0
    hspotcoff1(1)=2;
end
hspotcoff1(hf1)=hspotcoff1(hf1)+mm;
if hspotcoff1(hf1)>h
    hspotcoff1(hf1)=h;
end
vspotcoff1(1)=vspotcoff1(1)-mm;
if vspotcoff1(1)<0
    vspotcoff1(1)=2;
end
vspotcoff1(vf1)=vspotcoff1(vf1)+mm;
if vspotcoff1(vf1)>w
    vspotcoff1(vf1)=w;
end
disp(['drar grid line','coordinate y��',num2str(hf1),',coordinate x��',num2str(vf1)]);
%Drawgrid(p,h,w,hspotcoff1,vspotcoff1,hf1,vf1);
clear hspotcoff2;
clear vspotcoff2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%horizontal grid line optimizing
[hspotcoff2,hdmean,hdermean]=hgridadjust(hspotcoff1,0,0,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%vertial grid line optimizing
[vspotcoff2,vdmean,vdermean]=vgridadjust(vspotcoff1,0,0,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear hspotcoff1;
clear vspotcoff1;
clear hf2;
clear vf2;
et1=cputime-st;
disp(['time after optimizing',num2str(et1)]);
%darw the grid line
[s,hf2]=size(hspotcoff2);
[s,vf2]=size(vspotcoff2);
disp(['draw grid line: ','coordinate y��',num2str(hf2),',coordinate x��',num2str(vf2)]);
%Drawgrid(p,h,w,hspotcoff2,vspotcoff2,hf2,vf2);
%judge the gridding if correct or not, 1 is correct,0 is wrong
correct=iscorrect(hf,vf,hspotcoff2,vspotcoff2);
%output parameter
para=[et,et1,hf,vf,hf1,vf1,hf2,vf2,hdmean,vdmean,hdermean,vdermean,correct];
